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ABSTRACT 

We have developed a new three-dimensional general relativistic magneto- 
hydrodynamic (GRMHD) code, RAISHIN, using a conservative, high resolu- 
tion shock- capturing scheme. The numerical fluxes are calculated using the 
Harten, Lax, & van Leer (HLL) approximate Riemann solver scheme. The fiux- 
interpolated, constrained transport scheme is used to maintain a divergence-free 
magnetic field. In order to examine the numerical accuracy and the numer- 
ical efficiency, the code uses four different reconstruction methods: piecewise 
linear methods with Minmod and MC slope-limiter function, convex essentially 
non-oscillatory (CENO) method, and piecewise parabolic method (PPM) using 
multistep TVD Runge-Kutta time advance methods with second and third-order 
time accuracy. We describe code performance on an extensive set of test prob- 
lems in both special and general relativity. Our new GRMHD code has proven to 
be accurate in second order and has successfully passed with all tests performed, 
including highly relativistic and magnetized cases in both special and general 
relativity. 
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1. Introduction 

Both magnetic and gravitational fields play an important role in determining the evolu- 
tion of the matter in many astrophysical objects. In highly conducting plasma, the magnetic 
field can be amplified by gas contraction or shear motion. Even when the magnetic field is 
weak initially, the magnetic field grows short time scales and influences the gas dynamics of 
the system. This is particulary important for a compact object such as a black hole or a neu- 
tron star. Relativistic jets have been observed or postulated in various astrophysical objects, 
including active galactic nuclei (AGNs) (e.g.. Urry & Pavovani 1995; Ferrari 1998; Blandford 
2002), microquasars in our galaxy (e.g., Mirabel & Rodriguez 1999), and gamma-ray bursts 
(GRBs) (e.g., Zhang & Meszaros 2004; Piran 2005; Meszaros 2006). The most promising 
mechanisms for producing the relativistic jets involve the magnetohydrodynamic centrifugal 
acceleration and/or magnetic pressure driven acceleration from the accretion disk with the 
compact objects (e.g., Blandford & Payne 1982; Pukue 1990), or the extraction of rotating 
energy from the rotating black hole (Penrose 1969; Blandford & Znajek 1977). In addition, 
the differential rotation of the plasma in the disk raises the magnetorotational instability 
(MRI) , which plays an important role for the transportation of angular momentum from the 
disk due to the associated turbulence in accretion disks (Balbus & Hawley 1991, 1998). It 
is much more efficient than the internal disk transport. Magnetars are neutron stars with 
extremely large magnetic fields (~ 10^^ — lO^^G), as inferred from studies of anomalous X- 
ray pulsars and soft gamma-ray repeaters (e.g., Woods & Thompson 2006)). These intense 
magnetic fields are expected to affect the internal structure of the neutron star (e.g., Bocquet 
et al. 1995). The less intense fields of other neutron stars, ~ 10^^ — lO^^G, may also affect 
the internal structure. 

There has been much progresses in the analytical studies of relativistic astrophysical 
phenomena (e.g., Camenzind 1990; Fendt 1997; Fendt & Ouyed 2004). However, the problem 
is complex, involving time-dependent, three-dimensional dynamics of magnetized plasmas in 
the relativistic potential. Therefore analytical solutions are rather limited to simplified cases 
that are time stationary and spatially symmetric. Because of these, numerical experiments 
play an important role for complimenting theoretical works. 

A complete review of numerical approaches to relativistic hydrodynamics is given by 
Marti & Miiller (2003) and Font (2003). Numerical codes in special relativistic magnetohy- 
drodynamics (SRMHD) have been developed by a growing number of authors (van Putten 
1993; Koide et al. 1996; Komissarov 1999a; Balsara 2001; Koldoba et al. 2002; Del Zanna, 
Buccianti, & Londrillo 2003; Leismann et al. 2005; Mignone & Bodo 2006) and have been 
applied to the study of the relativistic jets (Koide et al. 1996, 1997; Koide 1997; Nishikawa 
et al. 1997, 1998; Kommisarov 1999b; Leismann et al. 2005; Mignone & Bodo 2006) and 
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pulsar wind nebulae (Kommisarov & Lyubarsky 2003, 2004; Bucciantini et al. 2003, 2004, 
2005, 2006; Del Zanna et al. 2004). In addition, Komissarov (1999a) and Balsara (2001) 
proposed a comprehensive set of tests to validate numerical codes in special relativity. The 
exact solutions of the Riemann problem in the SRMHD have been obtained by Giacomazzo 
& RezzoUa (2006). 

In order to investigate relativistic magnet orot at or s (RMRs), general relativistic mag- 
netohydrodynamics (GRMHD) codes with fixed spacctimes have been developed by some 
authors (Yokosawa 1993; Koide et al.l998; De Villiers & Hawley 2003; Gammie et al. 2003; 
Komissarov 2004; Anton et al. 2005; Anninos et al. 2005). These codes have been used 
to study the structure of accretion flows onto Kerr black holes and/or the formation of jets 
(Yokosawa 1995; Koide et al. 1998, 1999, 2000; Nishikawa et al. 2005; De Villiers et al. 
2003, 2005a; Hirose et al. 2004; Krolik et al. 2005; Hawley & Krolik 2006; McKinney & 
Gammie 2004; McKinney 2005, 2006), the Blandford-Znajek effect near the rotating black 
hole (Koide et al. 2002; Koide 2003; Komissarov 2005), and the formation of GRB jets in 
coUapsars (Mizuno et al. 2004a, 2004b; De Villiers et al. 2005b). 

Few attempts have been made to simulate relativistic MHD flows in dynamical space- 
times (e.g., Wilson 1975, Baumgarte & Shapiro 2003). Recently new codes capable of evolv- 
ing the Einstein-Maxwell-MHD equations have been developed by Duez et al. (2005), Shibata 
& Sekiguchi (2005) and Anderson ct al. (2006). These codes have been applied to investi- 
gate rotating neutron stars (Duez ct al. 2006a,b), and collapsing neutron stars for a central 
engine for short gamma-ray bursts (Shibata et al. 2006). These new GRMHD simulations 
with dynamical spacetimes will be applied to investigate various astrophysical systems such 
as merging neutron stars and black holes including gravitational waves and relativistic jets 
(e.g., Duez et al. 2005). 

Our previous GRMHD code developed by Koide has been applied to many high-energy 
astrophysical phenomena and showed pioneering results (Koide et al. 1998, 1999, 2000, 
2002; Koide 2003, 2004, 2006; Mizuno et al. 2004a,b; Nishikawa et al. 2005). However, the 
code cannot perform calculations in highly relativistic or highly magnetized regimes. The 
critical problem of the previous GRMHD code is that it cannot guarantee a divergence-free 
magnetic fleld. Even though we introduce a divergence- cleaning step in the code, it cannot 
perform long-term evolution, except in special cases. In order to overcome these numerical 
difficulties, we have developed a new, three-dimensional, GRMHD code called RAISHIN, for 
RelAtlviStic magnetoHydrodynamic sImulatioN. (RAISHIN is the ancient Japanese god of 
lightning.) It uses a conservative, high resolution shock-capturing scheme. 

The structure of the paper is as follows. In section 2 we describe the basic equations 
of GRMHD including the essentials of the 3-1-1 formalism, the description of the magnetic 
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field, the induction equation and the conservation equations of particle number, and the 
stress-energy tensor in conservative form. In section 3 we describe the basic algorithm of 
our new GRMHD code based on the recent modern techniques such as HLL approximate 
Riemann solver, various reconstruction methods and flux-interpolated constrained transport. 
We show the performance of the code on a series of test problems in section 4. The summary 
and conclusions are given in section 5. 

The new simulation results of jet formation using this code will be reported separately 
(Mizuno et al. 2006). 



2. Basic Equations 

2.1. Spacetimes and observers 

We investigate the evolution of a magnetized fluid in the background spacetime metric 
of a black hole written as 

ds'^ = -a^di^ + ^ij{dx' + p'dt) {dx^ + l3Ut), (1) 

where a, and are the lapse function, shift vector, and spatial metric respectively. A 
natural observer associated in the metric given by Eq. (1) is the so-called Eulerian observer 
with four- velocity perpendicular to the hypersurfaces of constant t at each event in the 
spacetime. The covariant and contravariant components of are given by 

71^ = 1(1,-/30, (2) 
a 

and 

n^ = (-a,0,0,0), (3) 

respectively. 

The comoving observer follows the fluid motion with four-velocity w''. The three- velocity 
of the fluid as measured by the Eulerian observer can be given as 

= (4) 

where —rinU^ = W is the relative Lorentz factor between and and /i^,^ = (7^,^ + n^rii, is 
the projector tensor onto the hypersurface orthogonal to n^. The spatial component of the 
projector tensor is hij — 7^^. Eq (4) is written as 
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and 

Vi = Ui/-f. (6) 
The Lorentz factor satisfies the following relation 

7 = — , — au^, (7) 



where v"^ — jij 



2.2. Evolution of the electromagnetic fields 

The electromagnetic field in general relativity is described by the Faraday electromag- 
netic tensor F'^'^. This tensor is related to the electric field E^^ and magnetic field B^^ 
measured by the Eulerian observer, 

F^'^ = nf'E'' - n'^Ef' + (8) 

where e^^""' — {—gy/'^[yijLua], g is the determinant of the four- metric {g — detgui) and e"^^"^ 
is the antisymmetric Levi-Civita symbol. Both electric and magnetic fields are orthogonal 
to n'^ {Ei^rifj, — B^Ufj, — 0), and can be expressed as 

= F^"'n^, (9) 

and 
where 

p*t^^ = -e^'^'^F^x (11) 

is dual of the electromagnetic field tensor. 

Wc adopt the ideal MHD condition and assume that the fluid is a perfect conductor. 
In this case the fluid has inflnite conductivity and in order to keep the current finite, the 
conduction current must vanish 

F^-^M^ = (12) 

which means that the electric field in the rest frame of the fiuid is zero. The electric and 
magnetic fields measured by a comoving observer are 



(13) 
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and 

B^" = u^F*''^'. (14) 

The ideal MHD condition (12) satisfies the condition that the electric field observed by a 
comoving observer becomes zero (£''' = 0). Bi^ is orthogonal to u^, i.e. li^-B'* = 0. The 
electromagnetic tensor can be expressed by the terms of B^^ as 

F^""" = u-ye'"""'B,. (15) 

The dual of the electromagnetic field tensor is also obtained by 

F*^"' = B^'u'' - B^'ui'. (16) 

Since B^^ is orthogonal to ti^, we have hi^B^ = B^. From Eqs. (10) and (16) 

h'^B'' = -nxu^B^. (17) 

Therefore 

^ _{V^. (18) 
The time and space components of Eq. (18) are given by 

B^ = (19) 
a. 

and 

B' = ^ ""f ■ (20) 

The evolution equation for the magnetic field can be obtained in conservation form from 
the dual of Maxwell's equation 

F^.,x + Fa^,. + = (21) 

in a coordinate basis, 

F*r = -^^{v^F*n-o. (22) 

Since ^/—g = 0:^/7 and F*** = B'^/a, the time component of Eq. (22) gives the divergence- 
free magnetic field constraint 

^ ^ i^B^) = 0, (23) 



and the spatial components of Eq. (22) gives the induction equation 

^ ^.X^B^) + ^^[^{u^B^-u^B^)]^0. (24) 
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It follows from Eq. (20) that 

where — — (3^ a. Therefore the induction equation can be written as 



g dx ' 



0. 



(25) 



(26) 



2.3. Conservation Equations 

The evolution equations for matter can be expressed as the local conservation laws 
for particle number and energy-momentum. The particle number conservation equation is 
written as 

{pu^;, = 0. (27) 
Here p is the rest-mass density. In a coordinate basis we rewrite this as 

1 d 



(V^pu^ = 0, 



and in 3-1-1 formalism 



-gDv') = 0, 



\f—9 dt \f~9 dx'^ 

where D — 7p. 

The energy-momentum conservation equation is given by 

where T^^ is the energy-momentum tensor. In a coordinate basis we rewrite this as 

1 d 



(28) 



(29) 



dx 



and in 3-1-1 formalism 

1 d 



-gdt 



1 d 

—= 7 

dx^ 



gT^ - r^^T- = 0, 



(30) 



(31) 



(32) 



where F^^, is the Christoffel symbol. The energy-momentum tensor for a system containing 
a perfect fluid and an electromagnetic field is the sum of a fluid part, 



TZa^phu>'u''+pg^\ 



(33) 
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where p is pressure, h is the specific enthalpy, defined hyh — 1 + u + p/p and u is internal 
energy and an electromagnetic part, 

T^M = ^ {f''Fx - l^FapF'^^^ . (34) 

In the ideal MHD condition, T^^j can be expressed simply in term of a magnetic 4-vector 
If = S/^/x/i^ as (e.g., Eq. (15)) 

= (35) 

where — b^b^,. Hence the energy-momentum tensor is given by 

Tf"" = {ph + b^)uV + (^P + ^) 9"" - b'^b''. (36) 

The spatial components of the energy-momentum conservation equation give the momentum 
equation 

where Si is the momentum density of the magnetized fluid 

Si = aTl = {ph + b'^h^Vj - ab%. (38) 

The time component of the energy-momentum conservation equation gives the energy equa- 
tion 

1 d 

' ^ ^ - - ' ' — — - — ' (39) 



where r is the total energy density 



r = a^r« -D^{ph + b')Y -{p + b'/2) - a\by - D. (40) 

To complete the system of equations, it remains to specify the equation of state (EOS). 
In this paper we adopt a F-law EOS 

p = (r - i)pu, (41) 

where V is adiabatic index. 

In summary, the evolution equations of GRMHD can be written in the following general 
form 

at <jx^ 
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where the quantities U, F, and S are 



D 

Si 

T 



(43) 



and 



where 0^ = (0, 0, 0)^ 



S = 



yiB3 - v^B' 



n^liu ( dgi,i per '\ 

0^ 



(44) 



(45) 



3. Numerical Scheme 

Numerical evolution of the GRMHD equations involves determining the fundamental 
MHD variables P = {p,p,v\ B^), called the "primitive" variables, at future times, given 
initial values of Pq. The evolution equations of GRMHD are written in conserved form 
(e.g., Eqs. (42)-(45)). They give the time derivatives for "conserved" variables U(P) = 
{D, Si, T, 5*) in terms of source variables S(P) and the divergence of flux variables F; 

There are several ways to evolve the GRMHD equations numerically. Conservative 
schemes evolve U with equation (46) at each time step. The advantage of these schemes 
is that highly accurate shock-capturing schemes can be apphed to the GRMHD equations. 
The disadvantage is that these schemes must recover P by numerically solving the system 
of equations, U = U(P), after each time step. This can be complicated. Conservative 
schemes for GRMHD have been developed by Koide et al. (1999), Koide (2003), Gammie et 
al. (2003), Komissarov (2004), Duez et al.(2005), Shibata & Sekiguchi (2005) and Anton et 
al.(2006). 
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On the other hand non-conservative schemes such as ZEUS (Stone & Norman 1992) 
evolve variables which are more simple with respect to P, but whose evolving equations are 
not the same form of Eq. (46). In such schemes, high resolution shock-capturing schemes 
cannot be applied and artificial viscosity must be used for handling discontinuities. The 
advantage of these schemes is they solve the internal energy equation instead of the energy 
equations (Eq. (39)). This is an advantage in regions where the internal energy is small 
compared to the total energy (such as highly supersonic flows). Moreover, the recovery of 
P is fairly straightforward. Non-conservative schemes for GRMHD have been developed by 
De Villiers & Hawley (2003) and Anninos et al. (2005). 

Our GRMHD code described in detail in the following section employs conservative 
schemes to solve the three-dimensional GRMHD equations on uniform and non-uniform 
grids in each spatial direction (method of lines). To maintain flexibihty our GRMHD code 
is programmed to allow for different boundary conditions, different coordinates (Cartesian, 
Cylindrical and Spherical in RMHD and Boyer-Lindquist coordinates in both non-rotating 
and rotating black holes), different spatial reconstruction algorithms, different time advance 
algorithms, and different recovery schemes. 



In order to improve the spatial accuracy of the code, we interpolate the primitive vari- 
ables within computational zones. These reconstructed variables are used to compute the 
fluxes F. For simplicity, we will consider the one-dimensional case. The generalization to 
three dimensions is straightforward. The primitive variables to the left and right of the grid 
cell interface arc P^ = Pj+i/2-e and Pr = Pj+i/2+e respectively. We have implemented 
several reconstruction schemes for computing Pl and Pr. 

1) Piecewise linear method (PLM) reconstruction 

We use the minmod slope-limited linear interpolation method and Monotonized central 
(MC) slope-hmited hnear interpolation method (van Leer 1977). These methods give 
results with the second-order accuracy at smooth regions and switch to flrst-order at 
local extrema. For given primitive variables a, values of a at the left and right of the 
grid cell interface ql and qr are computed according to 



3.1. The reconstruction step 



aL = ai + Va,;A,x72 



(47) 



qr = tti+i - Vai+iAa;/2 
Here, Va is the slope-hmited gradient of a. 



(48) 
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In the minmod slope limited linear interpolation method, 



Va = Ax ^minmod{ai^i — a^, 




(49) 



minmod{a, b) = 



sign{a)'min{\a\, \b\) 







if ah < 0, 
otherwise. 



(50) 



In the MC slope-limited linear interpolation method, 



MC{a, b) = 



Va = Ax ^MC{ai+i — a-i, ai — aj_i), 

r ifa6<0, 
\ sign{a)min{2\a\,2\b\,\a + b\/2) otherwise. 



(52) 



(51) 



2) Convex, essentially non-oscillatory (CENO) reconstruction 

In this scheme, polynomial (quadratic) interpolation is used to obtain the primitive 
variables of the left and right of the grid cell interface. In smooth regions, these values 
are accurate to third order in Ax. The scheme becomes first order at local extrema. 
The details of this scheme are written in Del Zanna & Bucciantini (2002) and Del 
Zanna et al. (2003). 

3) Piecewise parabohc method (PPM) reconstruction 

In this scheme, the quadratic polynomial interpolation is used to obtain the primitive 
variables to the left and right of the grid cell interface. These reconstructed values are 
then modified such that the parabohc profile defined by a^,, and is monotonic 
inside the grid cell. The modified interpolated values at the grid cell interfaces define 
local Riemann problems. In the regions near the contact discontinuities, the interpo- 
lation procedure is modified slightly to account sharp jumps. In the vicinity of the 
local extrema, the scheme switches to a piecewise constant approximation in order to 
avoid post shock oscillations. In smooth regions, these values are accurate to fourth 
order in Ax. The scheme becomes first order at local extrema. The details of this 
scheme are written in Colella & Woodward (1984) and the relativistic version of the 
PPM algorithm is written in Marti' and Miiller (1996). 



To calculate the numerical flux, we use the HLL (Harten, Lax, and van Leer) approx- 
imate Riemann solver (Harten et al. 1983). The HLL approximate Riemann solver is one 
of the simplest shock-capturing schemes because it does not require the eigenvector of the 



3.2. The Riemann solver step 
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characteristic matrix. However, when it couples with a high-order reconstruction scheme, 
it has been shown to perform with an accuracy comparable to more sophisticated solvers 
in shock tube problems (Del Zanna & Buccianitni 2002; Lucas-Serrano et al. 2004). To 
calculate the HLL fluxes, one only needs to know a maximum left-going wave speed c_|_ and 
a maximum right-going wave speed c_ at the both sides of the grid cell interface. 

Prom and P^,, we calculate the maximum left-going wave speeds c± the maximum 
right-going wave speeds c-t /j, the fluxes F/j = F(P^) and = F(Pi), and conserved 
variables IJr = IJ(Pr) and Ul = U(Pl). Defining c^ax = max{0,c+^R,c+^L) and Cmin = 
—min{0, C-^R, C-^l), the HLL flux is given by 

-r-i CminFiJ -|- CmaxF^ CminCmax(Ui^ /ro\ 

*i+i/2 = ■ (53) 

We calculate the wave speeds cj. by the same method as Gammie et al. (2003) and 
Duez et al. (2005). The HLL approximate Riemann solver requires only the maximum wave 
speeds in either direction along the three coordinate axes. To determine the wave speeds in 
the direction, one solves the dispersion relation for MHD waves with wave vectors of the 
form 

k^^(-u;,h,0,0). (54) 

The wave speed is the phase speed u/ki. The speeds along the x'^ and x^ direction are 
calculated in a similar way. The full dispersion relation for fast and slow modes is a fourth- 
order polynomial. It is difficult to solve the full dispersion relation. Gammie et al. (2003) 
has proposed replacing the full dispersion relation by the simpler approximate expression 
which overestimates the maximum speeds by a factor of < 2 (it makes the evolution more 
diffusive but more stable). In the frame comoving with the fluid, the approximate dispersion 
relation for slow and fast modes is written by 

u;'^K + cl{l-v\)]e, (55) 

where Cg = \/^vl (ph) is the sound speed, va is the Alfven speed, v\ = 6^/ {ph + h^) . 



3.3. Constrained transport 

High resolution shock-capturing scheme can successfully solve many problems involving 
various kinds of discontinuities. However these schemes do not guarantee V • B = in 
multidimensional simulations. Therefore we need to use constrained transport schemes to 
evolve the induction equation while maintaining V • B = 0. 
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Several approaches and schemes have been proposed to maintain V • B = (e.g., Toth 
2002). We use the flux- interpolated, constrained transport (flux-CT) scheme introduced 
by Toth (2002) and used by Gammie et al. (2003) and Duez et ah (2005). This scheme 
involves replacing the numerical flux of the induction equation computed at each point with 
linear combinations of the numerical fluxes computed at that point and neighboring points. 
The merit of this scheme is that it is naturally coupled with a Godunov-type scheme (Toth 
2002; Gammie et al 2003). It does not need the staggered variables which are needed in the 
constrained transport scheme by Evans & Hawley (1988), and the advantage of its usage 
with higher order implementation of the divergence free condition is discussed in Londrillo 
& Del Zanna (2004). 



In the time advance step, we get the updated values of the conserved variables at the 
next time levels (U"''"^). We use a multistep TVD Runge-Kutta (RK) method developed by 
Shu & Osher (1988) that can provide second (RK2) and third (RK3) order accuracy in time. 
The explicit form of the algorithms is: 

1. Prediction step (common for RK2 and RK3): 



3.4. Time advance step 



=U" + AiL(U"). 



(56) 



2. Depending on the accuracy of 



the time advance scheme do: 



RK2: 



U 



■n+l 



-[6U" + UW + AiL(uW)], 



(57) 



a 



in this case a — 2 and b — 1. 



RK3: 



-[6U" + UW + AtL(UW)], 



U(2) 



(58) 



a 



Un+l 




(59) 



in this case a — A and b — 3. 
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3.5. Recovery of primitive variables 

The conservative MHD scheme for GRMHD requires a method for transforming between 
conserved variables U and primitive variables P. The forward transformation (from primitive 
to conserved variables) has a closed-form solution, but the inverse transformation (from 
conserved to primitive variables) requires the solution of a set of five nonlinear equations. 
Having computed U at the new timestep, we must calculate primitive variables P at the 
new time. 

For the recovery we prepare two methods, Koide's 2D method (Koide et al. 1999) and 
Noble's 2D method (Noble et al. 2006). 



1. Koide's 2D method 

In Koide et al. (1999) they solve two nonlinear, simultaneous algebraic equations with 
two independent variables x = 7 — 1 and y = 7(v • B). In these equations they use 
r-law EOS (see Eq. (41)) 

r r 1' 

x{x + 2) TRx'^ + {2rR -d)x + rR-d + e + -y'^ 

= (Fx' + 2Vx + \f{f{x + 1)' + 2ay + 2axy + g^yX (60) 



V{R - g^)x^ + {2VR - 2r/ - d)x ^VR - d^ e - ^ -y^ 
= (T(x + l)(ra;2 + 2ra; + l), 



(61) 



where i? = D + r, rf=(r- \)D, e = (1 - r/2)Bl / = S, 5 = B, and a = B ■ S. 
Note that, in the absence of the magnetic field 5*, Eq (60) reduces to the well-known 
relativistic hydrodynamic one derived by Duncan and Hughes (1994), and Eq (61) 
becomes a trivial equation. These algebraic equations are solved at each grid cell 
using a 2-variable Newton-Raphson iteration method. The primitive variables then 
are calculated easily from x, y, S, r, and B, using 



(62) 

^ r.„;., , ' " , (63) 

+ {t+p + 57272 + (y/7) 72} ■ ^^^^ 
This method is identical to that used in special relativistic MHD simulations (Koide, 
Nishikawa, & Mutel 1996; Koide 1997). 



and 



7 = 1 + a:, 

(r - 1)[t -xD -(2- 1/7^)^72 + (t//7)72] 
Yix{x + 2) + 1] 

S + (Z//7)B 
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2. Noble's 2D method 

In Noble et al. (2006) they have tried six numerical methods for performing the inverse 
transformation and discuss the mathematical properties of them. We use one method 
which is recommended in Noble et al. (2006). This method solves two simple algebraic 
equations simultaneously for two independent variables W = 7^/1 and v'^ 

5=.(»y + B>^_(^!±£^)(^. (66) 

and ^ 

T^^{l + v') + ^ + W-D- p{u, p). (66) 

and 



Note that in Eq. (66) we can use any EOS. If we adopt a F-law EOS (Eq. (41)), Eq. (66) 
can be written as 

T^^{l + v') + ^ + W-D- (^[(1 - v')W - p]^ . (67) 

These algebraic equations arc solved at each grid cell using a 2-variable Newton-Raphson 
iteration method. From W and v"^ the primitive variables p and p (or u) are calculated easily. 



4. Code Tests 



In this section, we test the capabilities of our new GRMHD code. The tests are non- 
relativistic, special relativistic and general relativistic. 



4.1. Lineeir Alfven wave propagation 

The first test considers the propagation of a small- amphtude Alfven wave in one di- 
mension of Cartesian coordinates. The initial conditions are p — 1.0, p — 1.0, — 0.0, 
vy = Acos{kx), = 0.0, = So, By = -Bo{A/va) cos{kx), and = 0.0, where A; = 27r 
and A is the amphtude. We use the parameters Bq — 1.0 and A — 0.01. The fluid satisfies 
a F-law EOS with F = 4/3. The computational domain is < x < 1.0 and the boundary 
condition is periodic. 

The simulation runs for a single wave period 27r/a; {tend — Stt/o;), so that a perfect 
scheme would return to its initial state. We measure the Li norm of the difference between 
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the final state and the initial state for each primitive variables Q such as 

I'M 

Li{5Q) = / \Q{t = 0) - Q = Und)\dx (68) 
Ji=i 

as a function of the computational zone number N. 

Figure 1 shows the Li norm of the error in {Q = in Eq (68)) for runs using the MC 
slope limiter, minmod slope limiter, CENO, and PPM reconstructions as the computational 
zone number N is increased. Courant number is 0.5 in all simulations. All reconstruction 
schemes show that the global order of convergence of the code is second order in small 
computational zone number N and tend to flatter than second order convergence. Even 
we use the higher order reconstruction schemes such as CENO or PPM reconstructions, 
the global order of convergence of the code is second order. It is because we use the fiux- 
CT scheme to maintain divergence-free magnetic field. It has second-order accuracy. Even 
all reconstruction schemes show the second order convergence, MC slope-limiter and PPM 
reconstructions are more accurate than minmod slope limiter and CENO reconstruction 
schemes. 



4.2. Relativistic MHD shock-tube tests 

Shock-tube tests are the most basic test problems for MHD (HD) codes. Large sets 
of test-problems in relativistic MHD have been performed over the years (i.e., Komissarov 
1999a; Balsara 2001). In some test problems the exact solution of the Riemann problem in 
relativistic MHD has been calculated by Giacomazzo & RezzoUa (2006). 

We perform eight simulations of ^ cases which exact solutions obtained by Giaco- 
mazzo & Rezzolla (2006). Therefore we can compare the simulation results with the exact 
solutions directly. All tests start with discontinuous initial data at a; = (see Table 1) and 
with homogenous profiles on either side in Cartesian coordinates. We simulate from t = 
to t = tfinai with different reconstruction schemes, where tgnai is specified in Table 1 for each 
case. The fluid satisfies a F-law EOS. In all cases we use 400 computational zones with a 
Courant factor of 0.5. 

The results of the simulations arc shown in Fig 2 -17. All results show good agreement 
with the exact solutions. However, some small discontinuities and large shocks cannot be 
resolved exactly. This means that we need more computational zones to resolve all small 
discontinuities and large shocks exactly. Generally, the minmod slope hmiter and CENO re- 
constructions are more diffusive than the MC slope limiter and PPM reconstructions because 
of the properties of the minmod function. On the other hand, although the MC slope hmiter 
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and PPM reconstructions can resolve sharp discontinuities well, some small oscillations are 
seen at the discontinuities. The PPM reconstruction is the most accurate scheme to detect 
the discontinuities. 

Our previous GRMHD code of Koide (2003) could not handle some extreme cases of 
relativistic MHD shock-tube tests shown in Table 1 such as Kommissarov: Shock- Tube testl, 
Balsara Test2 and Balsara Tcst3 for large discontinuities of the pressure and magnetic field, 
Kommissarov: Collision Test and Balsara TestS for the highly relativistic flow even using 
different recovery methods such as Noble 2D method. However the new GRMHD code 
successfully handles the relativistic MHD shock-tube tests in Table 1. Therefore the new 
GRMHD code has substantial improvements over our previous GRMHD code (e.g., Koide 
2003) and can operate in a regime with large discontinuities of physical quantities (4 order 
difference of pressure in Komissarov: Collision Test and Balsara Test3), strong magnetic 
field {P < 0.004 and a > 570 in Balsara Test 3, where P — Pgas/Pmag, Pmag — ^'^/2, and 
a = l&^l/p) and highly relativistic flow (7 > 22 in Balsara Test 4). The hmitation to handle 
the regimes of high Lorentz factor and of highly magnetized depends on the schemes to solve 
the GRMHD equations. 



4.3. Magnetized Bondi flow 

Next we check the code to verify it numerically maintains the time-dependent system of 
equations describing the stationary, spherically symmetric accretion of a perfect fluid onto 
a Schwarzschild black hole in the presence of a radial magnetic field. Spherically symmetric 
accretion (Bondi flow) onto the fixed background of the Schwarzschild black hole has an 
analytic solution (e.g., Shapiro & Teukolsky 1983) that can be compared with the results of 
our code. This test has been used by several authors (De Villiers & Hawley 2003; Gammie 
et al. 2003; Duez et al. 2005; Shibata & Sekiguchi 2005; Anton et al. 2006) in the validation 
of their GRMHD codes. 

The initial setup consists of a perfect fluid which satisfles a F-law EOS with F = 4/3. 
The analytic solution is calculated in a manner similar to Koide et al. (1999) with constant 
parameter H — poho = 1.3 and po = 1-0. The critical point of free-falling flow is located 
at rs = 3.0. The radial magnetic held component is chosen to satisfy the divergence-free 
condition. Its strength is determined by the parameter Bq at the critical radius of the flow. 
These initial conditions are evolved in time in a uniform radial gird covering the region 
l.lrs < r < 20rs. 

Figure 18 shows the Li norm of the radial velocity {Q — in Eq. (68)) of the 
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difference between the exact solution and the final state {tend — SOts, Ts — rs/c ) in the 
Bq — 0.001 -^poc^ case. The calculation of Li norm is taken in the region 1.5rs < r < ISrs 
to exclude boundary effects. The MC slope limiter and minmod slope limiter reconstructions 
show that the global order of convergence of the code is second order. The MC slope limiter 
reconstruction is more accurate than the minmod slope limiter reconstruction. 

5. Summciry and Conclusions 

We have developed a new three-dimensional GRMHD code, RAISHIN, by using a con- 
servative high resolution shock-capturing scheme. The numerical fluxes are calculated using 
the HLL approximate Riemann solver scheme. The flux-interpolated, constrained transport 
scheme is used to maintain a divergence-free magnetic field. Several reconstruction and time 
advance schemes can be chosen for the numerical accuracy and computational resources. 

We have described several test problems in both special and general relativity. They 
have shown significant improvements over our previous GRMHD code (Koide 2003). Our 
new GRMHD code can perform in the regimes of high Lorentz factor (7 > 20) and high 
magnetic field (cr > 550), and in the presence of large discontinuity of density, pressure and 
magnetic field. We have compared the results of several reconstruction schemes. The code 
is second-order accurate even when we use the higher order reconstruction schemes such 
as CENO and PPM. Nevertheless, higher-order reconstruction schemes can provide more 
accurate results for some applications. The PPM reconstruction scheme allows the well- 
resolved detection of sharp discontinuities. We found the limitation to handle the regimes of 
high Lorentz factor and of highly magnetized depends on the schemes to solve the GRMHD 
equations. 

We have performed several simulations of non-rotating and rotating black hole systems 
with a thin accretion disk (Mizuno et al. 2006). The simulation results show the formation of 
jets driven by the Lorentz force and the gas pressure. It appears that the rotating black hole 
creates an additional, faster, and more collimated inner outfiow beside an outfiow generated 
by the rotating accretion disk in the non-rotating black hole. Thus, kinematic jet structure 
could be a sensitive function of the black hole rotation. 

Our new GRMHD code has proven to be accurate in second order and has successfully 
passed all applied tests including highly relativistic cases , and highly magnetized cases in 
both special and general relativity. We plan to apply this code to a number of high-energy 
astrophysical phenomena involving highly relativistic flows or compact objects with strong 
gravitational fields and magnetic fields. 
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Fig. 1. — Li norm of the error in for a linear Alfven wave propagation as a function of 
computational zone number, iV, for the MC slope limiter (plus), the minmod slope limiter 
(asterisk), convex CENO (open diamond) and PPM reconstructions (open triangle). The 
straight lines show the slope expected for second-order convergence (dashed hne) and third- 
order convergence (dash-dotted line). 
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Table 1. Relativistic MHD Shock-Tube Tests 



Test Type 




P 


P 










BV 




Komissarov: Shock Tube Testl left state 


1.0 


1000.0 


0.0 


0.0 


0.0 


1.0 


0.0 


0.0 


r = 4/3, tfinal = 1.0 


right state 


0.1 


1.0 


0.0 


0.0 


0.0 


1.0 


0.0 


0.0 


Komissarov: Collision Test 


left state 


1.0 


1.0 


5/^26 


0.0 


0.0 


10.0 


10.0 


0.0 


r = 4/3, tfinal = 1.2 


right state 


1.0 


1.0 


-5/^26 


0.0 


0.0 


10.0 


-10.0 


0.0 


Barsara Testl (Brio &: Wo) 


left state 


1.000 


1.0 


0.0 


0.0 


0.0 


0.5 


1.0 


0.0 


r = 2, tfinal = 0.4 


right state 0.125 


0.1 


0.0 


0.0 


0.0 


0.5 


-1.0 


0.0 


Barsara Test2 


left state 


1.0 


30.0 


0.0 


0.0 


0.0 


5.0 


6.0 


6.0 


r = 5/3, tfinal = 0.4 


right state 


1.0 


1.0 


0.0 


0.0 


0.0 


5.0 


0.7 


0.7 


Barsara Test 3 


left state 


1.0 


1000.0 


0.0 


0.0 


0.0 


10.0 


7.0 


7.0 


r = 5/3, tfinal = 0.4 


right state 


1.0 


0.1 


0.0 


0.0 


0.0 


10.0 


0.7 


0.7 


Barsara Test4 


left state 


1.0 


0.1 


0.999 


0.0 


0.0 


10.0 


7.0 


7.0 


r = 5/3, tfinal = 0.4 


right state 


1.0 


0.1 


-0.999 


0.0 


0.0 


10.0 


-7.0 


-7.0 


Barsara TestS 


left state 


1.08 


0.95 


0.40 


0.3 


0.2 


2.0 


0.3 


0.3 


r = 5/3, ifinal = 0.5 


right state 


1.00 


1.0 


-0.45 


-0.2 


0.2 


2.0 


-0.7 


0.5 


Generic Alfven Test 


left state 


1.0 


5.0 


0.0 


0.3 


0.4 


1.0 


6.0 


2.0 


r = 5/3, ifinal = 1-5 


right state 


0.9 


5.3 


0.0 


0.0 


0.0 


1.0 


5.0 


2.0 



Note. — 



Initial conditions for the relativistic shock-tube tests. 
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Fig. 2. — Simulation results of Komissarov shock-tube test 1 at time t — 1.0 using the MC 
slope hmiter (dotted lines) and the minmod slope hmiter (dashed hues) reconstructions. The 
solid lines arc the exact solution. The results are composed of a left-going fast rarefaction, 
a contact discontinuity, and a right-going fast shock. 
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Fig. 3. — Simulation results of Komissarov shock-tube test 1 at time t — l.Q using the CENO 
(dotted hues) and the PPM (dashed hues) reconstructions. The sohd hues are the exact 
solution. The results are composed of a left-going fast rarefaction, a contact discontinuity, 
and a right-going fast shock. 
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Fig. 4. — Simulation results of Komissarov collision test at time t — l.h using the MC slope 
limiter (dotted lines) and the minmod slope limiter (dashed hues) reconstructions. The 
solid lines are the exact solution. The results are composed of a left-going fast shock, of a 
left-going slow shock, a right-going slow shock, and a right-going fast shock. 
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Fig. 5. — Simulation results of Komissarov collision test at time t — 1.5 using the CENO 
(dotted lines) and the PPM (dashed hues) reconstructions. The sohd lines are the exact 
solution. The results arc composed of a left-going fast shock, of a left-going slow shock, a 
right-going slow shock, and a right-going fast shock. 
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Fig. 6. — Simulation results of Balsara test 1 at time t — 0.4 using the MC slope limiter 
(dotted lines) and the minmod slope hmiter (dashed hues) reconstructions. The sohd hues 
are the exact solutions. The results are composed of a left-going fast rarefaction, of a left- 
going slow shock, of a contact discontinuity, of a right-going slow shock and of a right-going 
fast rarefaction. 
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Fig. 7. — Simulation results of Balsara test 1 at time t — 0.4 using the CENO (dotted lines) 
and the PPM (dashed hues) reconstructions. The sohd lines are the exact solutions. The 
results are composed of a left-going fast rarefaction, of a left-going slow shock, a contact 
discontinuity, of a right-going slow shock, and a right-going fast rarefaction. 
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Fig. 8. — Simulation results of Balsara test 2 at time t — 0.4 using the MC slope limiter 
(dotted lines) and the minmod slope hmiter (dashed hues) reconstructions. The sohd hues 
are the exact solutions. The results are composed of two left-going fast and slow rarefactions, 
a contact discontinuity, and two right-going fast and slow shocks. 



-33- 




J 



I 



-nf fit u 



-«f «« «.< M 



-BhT u ^ 



(.g if^ e.* 



■!4a- ; \ — 



^1 I F.I 1:1 1 





^* A.f 



«.d e.i <L4 



-tlA -tf^ rAA-.- fl.t 



Fig. 9. — Simulation results of Balsara test 2 at time t — 0.4 using the CENO (dotted lines) 
and the PPM (dashed hues) reconstructions. The sohd lines are the exact solutions. The 
results are composed of two left-going fast and slow rarefactions, a contact discontinuity, and 
two right-going fast and slow shocks. 
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Fig. 10. — Simulation results of Balsara test 3 at time t — 0.4 using the MC slope limiter 
(dotted lines) and the minmod slope hmiter (dashed hues) reconstructions. The sohd hues 
are the exact solution. The results are composed of two left-going fast and slow rarefactions, 
a contact discontinuity, and two right-going fast and slow shocks. 
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Fig. 11. — Simulation results of Balsara test 3 at time t — 0.4 using the CENO (dotted lines) 
and the PPM (dashed lines) reconstructions. The solid hues are the exact solution. The 
results are composed of two left-going fast and slow rarefactions, a contact discontinuity, and 
two right-going fast and slow shocks. 
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Fig. 12. — Simulation results of Balsara test 4 at time t — 0.4 using the MC slope limiter 
(dotted lines) and the minmod slope hmiter (dashed hues) reconstructions. The sohd hues 
are the exact solution. The results are composed of two left-going fast and slow shocks, and 
two right-going fast and slow shocks. 
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Fig. 13. — Simulation results of Balsara test 4 at time t — 0.4 using the CENO (dotted 
lines) and the PPM (dashed hues) reconstructions. The solid lines are the exact solution. 
The results are composed of two left-going fast and slow shocks, and two right-going fast 
and slow shocks. 
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Fig. 14. — Simulation results of Balsara test 5 at time t — 0.5 using the MC slope limiter 
(dotted lines) and the minmod slope hmiter (dashed hues) reconstructions. The sohd hues 
arc the exact solution. The results are composed of a left-going fast shock, a left-going Alfven 
discontinuity, a left-going slow rarefaction, a contact discontinuity, a right-going slow shock, 
a right-going Alfven discontinuity, and a right-going fast shock. 
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Fig. 15. — Simulation results of Balsara test 5 at time t — 0.5 using the CENO (dotted 
lines) and the PPM (dashed hues) reconstructions. The solid lines are the exact solution. 
The results arc composed of a left-going fast shock, a left-going Alfvcn discontinuity, a left- 
going slow rarefaction, a contact discontinuity, a right-going slow shock, a right-going Alfven 
discontinuity, and a right-going fast shock. 
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Fig. 16. — Simulation results of Generic Alfven test at time t — 1.5 using the MC slope 
limiter (dotted lines) and the minmod slope hmiter (dashed lines) reconstructions. The solid 
lines arc the exact solution. The results are composed of a left-going fast rarefaction, a 
left-going Alfven discontinuity, a left-going slow shock, a contact discontinuity, a right-going 
slow shock, a right-going Alfven discontinuity, and a right-going fast shock. 
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Fig. 17. — Simulation results of Generic Alfven test at time t — 1.5 using the CENO (dotted 
lines) and the PPM (dashed hues) reconstructions. The solid lines are the exact solution. 
The results are composed of a left-going fast rarefaction, a left-going Alfven discontinuity, a 
left-going slow shock, a contact discontinuity, a right-going slow shock, a right-going Alfven 
discontinuity, and a right-going fast shock. 
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Fig. 18. — Li norm of the error in for a magnetized Bondi accretion flow as a function 
of computational zone number N for the MC slope limiter (plus) and the minmod slope 
limiter (asterisk) reconstructions. The straight hnes show the slope expected for second- 
order convergence (dashed line) and third-order convergence (dash-dotted line). 



